load obs

obs_only <- 
  read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
  mutate(subj = as_factor(subj),
         obs_degree = error,
         error = obs_degree * (pi/180))
## Parsed with column specification:
## cols(
##   subj = col_double(),
##   subj_index = col_double(),
##   stimulation = col_double(),
##   error = col_double()
## )
# obs_only <- sim_data %>%
#   unnest(subj_obs) %>%
#   select(subj_index = subj, 
#          error = obs_radian,
#          obs_degree, 
#          stimulation)

peek at data

pre
obs_only %>% 
  filter(stimulation == 0) %>%
  ggplot(aes(x = obs_degree)) +
  geom_histogram(binwidth = 10, aes(y=..density..)) + 
  geom_rug() + 
  geom_density(aes(y=..density..)) +  
  facet_wrap(vars(subj), ncol = 1)

post
obs_only %>% 
  filter(stimulation == 1) %>%
  ggplot(aes(x = obs_degree)) +
  geom_histogram(binwidth = 10, aes(y=..density..)) + 
  geom_rug() + 
  geom_density(aes(y=..density..)) +  
  facet_wrap(vars(subj), ncol = 1)

fit brms

source(glue("{params$model_dir_str}/model_prior.R"))

print(bprior_full)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.5)     b   intercept            circSD            
## 2   normal(0, 0.5)     b stimulation            circSD            
## 3   normal(0, 0.5)    sd   Intercept  subj      circSD            
## 4   normal(0, 0.5)    sd stimulation  subj      circSD            
## 5   normal(0, 1.5)     b   intercept             theta            
## 6   normal(0, 1.5)     b stimulation             theta            
## 7     normal(0, 1)    sd   Intercept  subj       theta            
## 8     normal(0, 1)    sd stimulation  subj       theta
iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

writeLines(
  make_stancode(bform_full, obs_only, family = vm_uniform_mix, prior = bprior_full, stanvars = stanvars),
  glue("{params$save_dir_str}/stan_code.txt")
)

model_fit <- brm(bform_full, obs_only, family = vm_uniform_mix, prior = bprior_full, stanvars = stanvars,
                 sample_prior = "yes",
                 warmup = warmup, iter = iter, cores = cores, chains = chains, 
                 control = list(adapt_delta = 0.99), inits = 0, 
                 file = glue("{params$save_dir_str}/obs_model_fit"))

print(model_fit)
##  Family: vm_uniform_mix 
##   Links: mu = identity; circSD = log; theta = logit; a = identity; b = identity 
## Formula: error ~ 0 
##          circSD ~ 0 + intercept + stimulation + (1 + stimulation || subj)
##          theta ~ 0 + intercept + stimulation + (1 + stimulation || subj)
##          a = -3.14
##          b = 3.14
##    Data: obs_only (Number of observations: 504) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~subj (Number of levels: 2) 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(circSD_Intercept)       0.36      0.25     0.02     0.96 1.00     3306
## sd(circSD_stimulation)     0.35      0.25     0.02     0.96 1.00     4334
## sd(theta_Intercept)        0.88      0.48     0.19     2.05 1.00     5352
## sd(theta_stimulation)      0.58      0.49     0.02     1.81 1.00     4306
##                        Tail_ESS
## sd(circSD_Intercept)       3035
## sd(circSD_stimulation)     3312
## sd(theta_Intercept)        4323
## sd(theta_stimulation)      3414
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## circSD_intercept       3.38      0.27     2.91     4.02 1.00     4535
## circSD_stimulation     0.26      0.29    -0.36     0.81 1.00     4995
## theta_intercept        0.10      0.62    -1.23     1.40 1.00     3782
## theta_stimulation     -0.11      0.59    -1.31     1.14 1.00     4571
##                    Tail_ESS
## circSD_intercept       4303
## circSD_stimulation     4531
## theta_intercept        3996
## theta_stimulation      3529
## 
## Family Specific Parameters: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## b_circSD_intercept       3.38      0.27     2.91     4.02 1.00     4535
## b_circSD_stimulation     0.26      0.29    -0.36     0.81 1.00     4995
## b_theta_intercept        0.10      0.62    -1.23     1.40 1.00     3782
## b_theta_stimulation     -0.11      0.59    -1.31     1.14 1.00     4571
##                      Tail_ESS
## b_circSD_intercept       4303
## b_circSD_stimulation     4531
## b_theta_intercept        3996
## b_theta_stimulation      3529
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

fit check

divergences

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

other

plot(model_fit, ask = FALSE)

* Notes

no divergences,

good ESS,

good rhat,

good fuzzy plots

plot posteriors

arrange samples

# compute summaries for plot

group_level_samples <- 
  spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
  mutate(
         # group level parameters
         circSD_pre_mean  = exp(b_circSD_intercept),
         circSD_post_mean = exp(b_circSD_intercept + b_circSD_stimulation),
         circSD_ES_mean   = circSD_post_mean - circSD_pre_mean,
         pMem_pre_mean    = inv_logit(b_theta_intercept),
         pMem_post_mean   = inv_logit(b_theta_intercept + b_theta_stimulation),
         pMem_ES_mean     = pMem_post_mean - pMem_pre_mean,
         # predicitve dist for group level parameters
         circSD_pre_pred  = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept)),
         circSD_post_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept) + 
                                rnorm(n(), b_circSD_stimulation, sd_subj__circSD_stimulation)),
         circSD_ES_pred   = circSD_post_pred - circSD_pre_pred,
         pMem_pre_pred    = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept)),
         pMem_post_pred   = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept) +
                                    rnorm(n(), b_theta_stimulation, sd_subj__theta_stimulation)),
         pMem_ES_pred     = pMem_post_pred - pMem_pre_pred
         ) %>% 
  select(-contains("b_"), -contains("sd_subj")) %>%
  pivot_longer(-contains("."), names_to = c("param", "stat"), names_pattern = "(.*)_(.*)", values_to = "value") %>%
  pivot_wider(names_from = stat, values_from = value)
  
group_level_summary <- 
  group_level_samples %>%
  group_by(param) %>%
  median_qi(.width = c(.5, .8, .95))



circSD_subj_samples <- 
  model_fit %>%
  spread_draws(b_circSD_intercept, b_circSD_stimulation, r_subj__circSD[subj, term]) %>%
  ungroup() %>%
  pivot_wider(names_from = term, values_from = r_subj__circSD, names_prefix = "offset_") %>%
  mutate(subj = subj,
        circSD_pre = exp(b_circSD_intercept + offset_Intercept),
        circSD_post = exp(b_circSD_intercept + offset_Intercept + b_circSD_stimulation + offset_stimulation),
        circSD_ES = circSD_post - circSD_pre) %>%
  select(-c(b_circSD_intercept, offset_Intercept, b_circSD_stimulation, offset_stimulation)) %>%
  pivot_longer(contains("circSD"), names_to = "param", values_to = "value") 

circSD_subj_summary <- 
  circSD_subj_samples %>%
  group_by(subj, param) %>%
  median_qi(.width = c(.90, .95))
## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed

## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed
pMem_subj_samples <- 
  model_fit %>%
  spread_draws(b_theta_intercept, b_theta_stimulation, r_subj__theta[subj, term]) %>%
  ungroup() %>%
  pivot_wider(names_from = term, values_from = r_subj__theta, names_prefix = "offset_") %>%
  mutate(subj = subj,
            pMem_pre = inv_logit(b_theta_intercept + offset_Intercept),
            pMem_post = inv_logit(b_theta_intercept + offset_Intercept + b_theta_stimulation + offset_stimulation),
            pMem_ES = pMem_post - pMem_pre) %>%
  select(-c(b_theta_intercept, offset_Intercept, b_theta_stimulation, offset_stimulation)) %>%
  pivot_longer(contains("pMem"), names_to = "param", values_to = "value") 

pMem_subj_summary <- 
  pMem_subj_samples %>%
  group_by(subj, param) %>%
  median_qi(.width = c(.90, .95))
## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed

## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed
group_level_samples %>%
  select(-pred) %>%
  group_by(param) %>%
  median_qi(.width = c(.95))
## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed
## # A tibble: 6 x 7
##   param           mean  .lower .upper .width .point .interval
##   <chr>          <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 circSD_ES     8.88    -9.26  38.7     0.95 median qi       
## 2 circSD_post  37.6     18.5   83.7     0.95 median qi       
## 3 circSD_pre   28.7     18.3   55.7     0.95 median qi       
## 4 pMem_ES      -0.0282  -0.283  0.239   0.95 median qi       
## 5 pMem_post     0.494    0.152  0.849   0.95 median qi       
## 6 pMem_pre      0.523    0.226  0.802   0.95 median qi

group level posteriors

circSD_p1 <- group_level_samples %>% 
  filter(str_detect(param, "circSD")) %>%
  ggplot() + 
  # posterior dist + interval for group mean
  geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) + 
  # posterior predictive distribution for group means
  stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
  # posterior medians for each parameter estimate per subj
  geom_point(data = circSD_subj_summary, aes(y = param, x = value), size = 2) +
  # decorations
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  coord_cartesian(xlim=c(-50, 150)) +
  labs(subtitle = "circSD: group level mean posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects", 
       x = "circSD", 
       color = "interval")

# group level pMem pre, post and ES plot

pMem_p1 <- group_level_samples %>% 
  filter(str_detect(param, "pMem")) %>%
  ggplot() + 
  # posterior dist + interval for group mean
  geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) + 
  # posterior predictive distribution for group means
  stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
  # posterior medians for each parameter estimate per subj
  geom_point(data = pMem_subj_summary, aes(y = param, x = value), size = 2) +
  # decorations
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  #coord_cartesian(xlim=c(-50, 150)) +
  labs(subtitle = "pMem: group level mean posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects", 
       x = "pMem", 
       color = "interval")

plot_grid(circSD_p1, pMem_p1, align = "hv", ncol = 1)

circSD pre, post, ES

# circSD pre: group level posteriors and subject posteriors

circSD_p2 <- 
  ggplot() + 
  # plot group mean circSD_pre posterior and subject circSD_pre posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "circSD_pre")) %>%
                              select(-pred, value = mean)
                            ,
                              circSD_subj_samples %>%
                              filter(str_detect(param, "circSD_pre")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot pre condition group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_pre")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject circSD_pre posterior medians in the prediction band
  geom_point(data = circSD_subj_summary %>% filter(param == "circSD_pre"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y = -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "circSD_pre: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
       x = "circSD",
       color = "interval")



circSD_p3 <-   
  ggplot() + 
  # plot group mean circSD_post posterior and subject circSD_post posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "circSD_post")) %>%
                              select(-pred, value = mean)
                            ,
                              circSD_subj_samples %>%
                              filter(str_detect(param, "circSD_post")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot post condition group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_post")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject circSD_post posterior medians in the prediction band
  geom_point(data = circSD_subj_summary %>% filter(param == "circSD_post"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y =  -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "circSD_post: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
       x = "circSD",
       color = "interval")


circSD_p4 <- 
  ggplot() + 
  # plot group mean circSD_ES posterior and subject circSD_ES posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "circSD_ES")) %>%
                              select(-pred, value = mean)
                            ,
                              circSD_subj_samples %>%
                              filter(str_detect(param, "circSD_ES")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot ES  group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_ES")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject circSD_ES posterior medians in the prediction band
  geom_point(data = circSD_subj_summary %>% filter(param == "circSD_ES"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y =  -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "circSD_ES: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \nES predictive dist of subjects",
       x = "detla circSD",
       color = "interval")


plot_grid(circSD_p2, circSD_p3, circSD_p4, ncol = 1, align = "hv")

pMem pre, post, ES

# pMem_ pre: group level posteriors and subject posteriors

pMem_p2 <- 
  ggplot() + 
  # plot group mean pMem__pre posterior and subject cpMem__pre posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "pMem_pre")) %>%
                              select(-pred, value = mean)
                            ,
                              pMem_subj_samples %>%
                              filter(str_detect(param, "pMem_pre")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot pre condition group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_pre")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject pMem_pre posterior medians in the prediction band
  geom_point(data = pMem_subj_summary %>% filter(param == "pMem_pre"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y = -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "pMem_pre: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
       x = "pMem",
       color = "interval")



pMem_p3 <-   
  ggplot() + 
  # plot group mean pMem_post posterior and subject pMem_post posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "pMem_post")) %>%
                              select(-pred, value = mean)
                            ,
                              pMem_subj_samples %>%
                              filter(str_detect(param, "pMem_post")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot post condition group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_post")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject pMem_post posterior medians in the prediction band
  geom_point(data = pMem_subj_summary %>% filter(param == "pMem_post"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y =  -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "pMem_post: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
       x = "pMem",
       color = "interval")


pMem_p4 <- 
  ggplot() + 
  # plot group mean pMem_ES posterior and subject pMem_ES posterior
  geom_halfeyeh(data = rbind(
                              group_level_samples %>% 
                              filter(str_detect(param, "pMem_ES")) %>%
                              select(-pred, value = mean)
                            ,
                              pMem_subj_samples %>%
                              filter(str_detect(param, "pMem_ES")) %>%
                              unite(param, param, subj) )
                , aes(y = param, x = value), .width = c(.90, .95)) + 
  # plot ES  group predictive distribution
  stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_ES")),
                 aes(y = param , x = pred),
                 .width = c(.5, .8, .95),
                 position = position_nudge(y = -0.15)
                   ) +
  # show modeled subject pMem_ES posterior medians in the prediction band
  geom_point(data = pMem_subj_summary %>% filter(param == "pMem_ES"),
             aes(y = param, x = value), 
             size = 2, 
             position = position_nudge(y =  -0.15)) + 
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  labs(subtitle = "pMem_ES: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \nES predictive dist of subjects",
       x = "detla pMem",
       color = "interval")
plot_grid(pMem_p2, pMem_p3, pMem_p4, ncol = 1, align = "hv")

group level joint posteriors

group_level_samples %>%
  pivot_wider(id_cols = contains("."), names_from = param, values_from = mean) %>%
  select(-contains(".")) %>%
  mcmc_pairs(off_diag_fun = "hex")
## Warning: Only one chain in 'x'. This plot is more useful with multiple
## chains.

plot posteriors w/ priors

# compute summaries for plot with priors

group_level_samples <- 
  spread_draws(model_fit, `(prior_)?(b|sd)_.*`, regex = TRUE) %>%
  mutate(
         # prior
         # group level parameters
         prior_circSD_pre_mean  = exp(prior_b_circSD_intercept),
         prior_circSD_post_mean = exp(prior_b_circSD_intercept + prior_b_circSD_stimulation),
         prior_circSD_ES_mean   = prior_circSD_post_mean - prior_circSD_pre_mean,
         prior_pMem_pre_mean    = inv_logit(prior_b_theta_intercept),
         prior_pMem_post_mean   = inv_logit(prior_b_theta_intercept + prior_b_theta_stimulation),
         prior_pMem_ES_mean     = prior_pMem_post_mean - prior_pMem_pre_mean,   
         # predicitve dist for group level parameters
         prior_circSD_pre_pred  = exp(rnorm(n(), prior_b_circSD_intercept, prior_sd_subj__circSD_Intercept)),
         prior_circSD_post_pred = exp(rnorm(n(), prior_b_circSD_intercept, prior_sd_subj__circSD_Intercept) + 
                                rnorm(n(), prior_b_circSD_stimulation, prior_sd_subj__circSD_stimulation)),
         prior_circSD_ES_pred   = prior_circSD_post_pred - prior_circSD_pre_pred,
         prior_pMem_pre_pred    = inv_logit(rnorm(n(), prior_b_theta_intercept, prior_sd_subj__theta_Intercept)),
         prior_pMem_post_pred   = inv_logit(rnorm(n(), prior_b_theta_intercept, prior_sd_subj__theta_Intercept) +
                                    rnorm(n(), prior_b_theta_stimulation, prior_sd_subj__theta_stimulation)),
         prior_pMem_ES_pred     = prior_pMem_post_pred - prior_pMem_pre_pred,
         
         # posteriors
         # group level parameters
         circSD_pre_mean  = exp(b_circSD_intercept),
         circSD_post_mean = exp(b_circSD_intercept + b_circSD_stimulation),
         circSD_ES_mean   = circSD_post_mean - circSD_pre_mean,
         pMem_pre_mean    = inv_logit(b_theta_intercept),
         pMem_post_mean   = inv_logit(b_theta_intercept + b_theta_stimulation),
         pMem_ES_mean     = pMem_post_mean - pMem_pre_mean,
         # predicitve dist for group level parameters
         circSD_pre_pred  = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept)),
         circSD_post_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept) + 
                                rnorm(n(), b_circSD_stimulation, sd_subj__circSD_stimulation)),
         circSD_ES_pred   = circSD_post_pred - circSD_pre_pred,
         pMem_pre_pred    = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept)),
         pMem_post_pred   = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept) +
                                    rnorm(n(), b_theta_stimulation, sd_subj__theta_stimulation)),
         pMem_ES_pred     = pMem_post_pred - pMem_pre_pred
         ) %>% 
  select(-contains("b_"), -contains("sd_subj")) %>%
  pivot_longer(-contains("."), names_to = c( "param", "stat"), names_pattern = "(.*)_(.*)", values_to = "value") %>%
  pivot_wider(names_from = stat, values_from = value)

group level posteriors

circSD_p1_wPrior <- 
  group_level_samples %>% 
  filter(str_detect(param, "circSD")) %>% 
  mutate(param = fct_relevel(param, c("prior_circSD_ES", "circSD_ES", "prior_circSD_pre", "circSD_pre", "prior_circSD_post", "circSD_post"))) %>%
  ggplot() + 
  
  # posterior dist + interval for group mean
  geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) + 
  # posterior predictive distribution for group means
  stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
  # posterior medians for each parameter estimate per subj
  geom_point(data = circSD_subj_summary, aes(y = param, x = value), size = 2) +
  
  # prior dist + interval for group mean
  #geom_halfeyeh(data = group_level_prior_param_samples, aes(y = param, x = values), .width = c(.90, .95)) +
  
  # decorations
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  coord_cartesian(xlim=c(-50, 150)) +
  labs(subtitle = "circSD: group level mean prior + posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects", 
       x = "circSD", 
       color = "interval")



# group level pMem pre, post and ES plot

pMem_p1_wPrior <- 
  group_level_samples %>% 
  filter(str_detect(param, "pMem")) %>%
  mutate(param = fct_relevel(param, c("prior_pMem_ES", "pMem_ES", "prior_pMem_pre", "pMem_pre", "prior_pMem_post", "pMem_post"))) %>%
  ggplot() + 
  
  # posterior dist + interval for group mean
  geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) + 
  # posterior predictive distribution for group means
  stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
  # posterior medians for each parameter estimate per subj
  geom_point(data = pMem_subj_summary, aes(y = param, x = value), size = 2) +
  
  # decorations
  scale_color_brewer() + 
  scale_x_continuous(breaks = pretty_breaks(10)) + 
  #coord_cartesian(xlim=c(-50, 150)) +
  labs(subtitle = "pMem: group level mean prior + posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects", 
       x = "pMem", 
       color = "interval")

plot_grid(circSD_p1_wPrior, pMem_p1_wPrior, align = "hv", ncol = 1)

* Notes

Main conclusion is that both group mean and subject predicitive dist show reductions in uncertainty, though lots of probability still over zero for both pMem and circSD ES.

There is a sign of a circSD increase, more data will make that clearer.

posterior predictive plot of errors

sim

conditions <- c(0,1)
nobs_per_cond_sim <- 1500

post_pred_fpath <- glue("{params$save_dir_str}/post_pred_obs.rds")

if (file.exists(post_pred_fpath)){
  
  post_pred_sim <- readRDS(post_pred_fpath)

}else{

  print("simulating")
  
  posterior_arranged <- 
    spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
    sample_n(2e3) %>%
    mutate(sim_num = 1:n(),
           nsubj = 1,
           nobs_per_cond = nobs_per_cond_sim) %>%
    select(-contains(".")) %>%
    select(sim_num,
           alpha0_mu = b_circSD_intercept,
           alpha0_sigma = sd_subj__circSD_Intercept,
           alphaD_mu = b_circSD_stimulation,
           alphaD_sigma = sd_subj__circSD_stimulation,
           beta0_mu = b_theta_intercept,
           beta0_sigma = sd_subj__theta_Intercept,
           betaD_mu = b_theta_stimulation,
           betaD_sigma = sd_subj__theta_stimulation,
           nsubj,
           nobs_per_cond) 
  
  post_pred_sim <- 
    posterior_arranged %>%
    mutate(
      # use draw_subj to sample nsubj_sim per sim using group-level parameter draws
      dataset = pmap(posterior_arranged, draw_subj),
      stimulation = list(stimulation = rep(conditions, each = nobs_per_cond_sim))) %>%
    
    # first unnest dataset, expanding by nsubj_sim and copying stimulation list to each subj
    unnest(dataset) %>%
    
    # then unnest stimulation, expanding by nobs_per_cond_sim*2
    unnest(stimulation) %>%
    
    # now use likelihood to simulation observations
    mutate(
      # evaluate and delink linear model on pMem
      pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
      
      # evaluate and delink linear model on circSD/kappa
      k = sd2k_vec(
        pracma::deg2rad(
          exp(subj_alpha0 + (subj_alphaD * stimulation)))),
      
      # use pMem to draw a 1 or 0 for each trial
      memFlip = rbernoulli(n(), pMem),
      
      # use k to draw from vonMises for each trial
      vm_draw = rvonmises_vec(1, pi, k) - pi,
      
      # draw from unif for each trial
      unif_draw = runif(n(), -pi, pi),
      
      # assign either vm_draw or unif_draw to each trial, depending on memFlip
      obs_radian = memFlip * vm_draw + (1 - memFlip) * unif_draw,
      
      # convert to degrees
      obs_degree = obs_radian * (180/pi)
    ) %>%
    select(-c(pMem, k, memFlip, vm_draw, unif_draw)) %>%
    nest(subj_obs = c(stimulation, obs_degree, obs_radian)) %>%
    nest(dataset = c(subj, subj_alpha0, subj_alphaD, subj_beta0, subj_betaD, nobs_per_condition, subj_obs))

  
  saveRDS(post_pred_sim, post_pred_fpath)
  
}
conditions <- c(0,1)
nobs_per_cond_sim <- 1500

post_predmean_fpath <- glue("{params$save_dir_str}/post_predmean_obs.rds")

if (file.exists(post_predmean_fpath)){
  
  post_predmean_sim <- readRDS(post_predmean_fpath)

}else{

  posterior_arranged <- 
    spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
    sample_n(2e3) %>%
    mutate(sim_num = 1:n(),
           nobs_per_cond = nobs_per_cond_sim ) %>%
    select(-contains(".")) %>%
    select(sim_num,
           subj_alpha0 = b_circSD_intercept,
           subj_alphaD = b_circSD_stimulation,
           subj_beta0 = b_theta_intercept,
           subj_betaD = b_theta_stimulation,
           nobs_per_cond) 

  post_predmean_sim <- 
    posterior_arranged %>%
    
    # add stimulation covariates
    mutate(
      stimulation = list(stimulation = rep(conditions, each = nobs_per_cond_sim))
    ) %>%
    
    # then unnest stimulation, expanding by nobs_per_cond_sim*2
    unnest(stimulation) %>%
    
    # now use likelihood to simulation observations
    mutate(
      
      # evaluate and delink linear model on pMem
      pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
      
      # evaluate and delink linear model on circSD/kappa
      k = sd2k_vec(
        pracma::deg2rad(
          exp(subj_alpha0 + (subj_alphaD * stimulation)))),
      
      # use pMem to draw a 1 or 0 for each trial
      memFlip = rbernoulli(n(), pMem),
      
      # use k to draw from vonMises for each trial
      vm_draw = rvonmises_vec(1, pi, k) - pi,
      
      # draw from unif for each trial
      unif_draw = runif(n(), -pi, pi),
      
      # assign either vm_draw or unif_draw to each trial, depending on memFlip
      obs_radian = memFlip * vm_draw + (1 - memFlip) * unif_draw,
      
      # convert to degrees
      obs_degree = obs_radian * (180/pi)
      
    ) %>%
    select(-c(pMem, k, memFlip, vm_draw, unif_draw)) %>%
    nest(dataset = c(stimulation, obs_degree, obs_radian)) 

  
  saveRDS(post_predmean_sim, post_predmean_fpath)
  
}

posterior data prediction, w/ group-level variance

#######################################################
# calculate histogram quantile mats from each condition

sim_subj_obs_hist_count <- function(dataset, condition = 0){
  
  dataset_obs <- dataset %>% 
    sample_n(1) %>%
    unnest(subj_obs) %>%
    filter(stimulation == condition) %>%
    select(obs_degree)
  
  breaks <- seq(-180, 180, 5)
  
  bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
  
  bincount_names <- glue("c{breaks[-1]}")
  
  names(bincount) <- bincount_names
  bincount_df <- data.frame(as.list(bincount))

  return(bincount_df)
  
}

make_quantmat <- function(sim_datasets, condition = 0){

  bincounts <- sim_datasets %>% 
  select(dataset) %>% 
  mutate(subj_hist_counts = map(dataset, sim_subj_obs_hist_count, condition)) %>% 
  select(-dataset) %>% 
  unnest(subj_hist_counts) %>%
  as_tibble()


  xvals <- seq(-177.5, 177.5, 5)
  probs <- seq(0.1,0.9,0.1)

  quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
  names(quantmat) <- paste0("p",probs)

  quantmat <- cbind(quantmat, xvals)

  for (iQuant in 1:length(probs)){
   quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
  }
    
  return(quantmat)
}

quantmat_cond0 <- make_quantmat(post_pred_sim, 0)
quantmat_cond1 <- make_quantmat(post_pred_sim, 1)


#######################################################
# calculate ecdf quantile mats from each condition

# unnest post_pred_sim, using only 1 subj/dataset
unnested <- 
  post_pred_sim %>%
  unnest(dataset) %>%
  group_by(sim_num) %>%
  sample_n(1) %>%
  ungroup() %>%
  unnest(subj_obs)

# calc quantiles mat for pre condition
ecdf_res_stim0 <- 
  unnested %>% 
  filter(stimulation == 0) %>%
  group_by(sim_num) %>% 
  group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))

stim0_ecdf_quantiles <- bind_cols(
  tibble(x_val = seq(-180, 180, 1)), 
  as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim0), probs = c(0.95, 0.5, 0.05 )))
  )

# calc quantiles mat for post condition
ecdf_res_stim1 <- 
  unnested %>% 
  filter(stimulation == 1) %>%
  group_by(sim_num) %>% 
  group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))

stim1_ecdf_quantiles <- bind_cols(
  tibble(x_val = seq(-180, 180, 1)), 
  as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim1), probs = c(0.95, 0.5, 0.05 )))
  )
# blues
b_light <- "#8C9BC4"
b_light_highlight <- "#A0ADCE"
b_mid   <- "#546BA9"
b_mid_highlight   <- "#7385B8"
b_dark  <- "#002381"
b_dark_highlight  <- "#2E4B97"

#betancourt reds
r_light <- "#DCBCBC"
r_light_highlight <- "#C79999"
r_mid   <- "#B97C7C"
r_mid_highlight   <- "#A25050"
r_dark  <- "#8F2727"
r_dark_highlight  <- "#7C0000"


#######################################################
# plot histogram(density) per condition

ggplot(quantmat_cond0, aes(x = xvals)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = r_light, alpha = 0.4) + 
  geom_line(aes(y = p0.5), color = r_dark, size = 1) + 
  geom_ribbon(data = quantmat_cond1, aes(ymax = p0.9, ymin = p0.1), fill = b_light, alpha = 0.4) + 
  geom_line(data = quantmat_cond1, aes(y = p0.5), color = b_dark, size = 1) + 
  scale_x_continuous(breaks=pretty_breaks(10)) + 
  labs(x = "error (degrees) [red = pre, blue = post]", 
       y = "count +/- quantile", 
       subtitle = glue("per-condition posterior pred dist (median line, 90% interval over {nrow(post_pred_sim)} sim datasets)\n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean + sd draw)")) + 
  theme_cowplot()

#######################################################
# plot ecdf per condition

ggplot() +
  geom_ribbon(data = stim0_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "red", alpha = 0.3) + 
  geom_line(data = stim0_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "red", size = 1) +
  geom_ribbon(data = stim1_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "blue", alpha = 0.3) + 
  geom_line(data = stim1_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "blue", size = 1) + 
  scale_x_continuous(breaks=pretty_breaks(10)) + 
  geom_hline(yintercept = seq(0, 1, 0.25), linetype = "dashed", alpha = 0.2) +
  labs(x = "error (degrees) [red = pre, blue = post]", 
       y = "cumulative prob.", 
       subtitle = glue("per-condition posterior pred cdf (median line, 90% interval over {nrow(post_pred_sim)} sim datasets) \n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean + sd draw")) + 
  theme_cowplot()

posterior data prediction, w/o group-level variance

#######################################################
# calculate histogram quantile mats from each condition

sim_single_subj_obs_hist_count <- function(dataset, condition = 0){
  
  dataset_obs <- dataset %>% 
    filter(stimulation == condition) %>%
    select(obs_degree)
  
  breaks <- seq(-180, 180, 5)
  
  bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
  
  bincount_names <- glue("c{breaks[-1]}")
  
  names(bincount) <- bincount_names
  bincount_df <- data.frame(as.list(bincount))

  return(bincount_df)
  
}

make_quantmat <- function(sim_datasets, condition = 0){

  bincounts <- sim_datasets %>% 
  select(dataset) %>% 
  mutate(subj_hist_counts = map(dataset, sim_single_subj_obs_hist_count, condition)) %>% 
  select(-dataset) %>% 
  unnest(subj_hist_counts) %>%
  as_tibble()


  xvals <- seq(-177.5, 177.5, 5)
  probs <- seq(0.1,0.9,0.1)

  quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
  names(quantmat) <- paste0("p",probs)

  quantmat <- cbind(quantmat, xvals)

  for (iQuant in 1:length(probs)){
   quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
  }
    
  return(quantmat)
}

quantmat_cond0 <- make_quantmat(post_predmean_sim, 0)
quantmat_cond1 <- make_quantmat(post_predmean_sim, 1)


#######################################################
# calculate ecdf quantile mats from each condition

# unnest post_pred_sim, using only 1 subj/dataset
unnested <- 
  post_predmean_sim %>%
  unnest(dataset)

# calc quantiles mat for pre condition
ecdf_res_stim0 <- 
  unnested %>% 
  filter(stimulation == 0) %>%
  group_by(sim_num) %>% 
  group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))

stim0_ecdf_quantiles <- bind_cols(
  tibble(x_val = seq(-180, 180, 1)), 
  as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim0), probs = c(0.95, 0.5, 0.05 )))
  )

# calc quantiles mat for post condition
ecdf_res_stim1 <- 
  unnested %>% 
  filter(stimulation == 1) %>%
  group_by(sim_num) %>% 
  group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))

stim1_ecdf_quantiles <- bind_cols(
  tibble(x_val = seq(-180, 180, 1)), 
  as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim1), probs = c(0.95, 0.5, 0.05 )))
  )
# blues
b_light <- "#8C9BC4"
b_light_highlight <- "#A0ADCE"
b_mid   <- "#546BA9"
b_mid_highlight   <- "#7385B8"
b_dark  <- "#002381"
b_dark_highlight  <- "#2E4B97"

#betancourt reds
r_light <- "#DCBCBC"
r_light_highlight <- "#C79999"
r_mid   <- "#B97C7C"
r_mid_highlight   <- "#A25050"
r_dark  <- "#8F2727"
r_dark_highlight  <- "#7C0000"


#######################################################
# plot histogram(density) per condition

ggplot(quantmat_cond0, aes(x = xvals)) + 
  geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = r_light, alpha = 0.4) + 
  geom_line(aes(y = p0.5), color = r_dark, size = 1) + 
  geom_ribbon(data = quantmat_cond1, aes(ymax = p0.9, ymin = p0.1), fill = b_light, alpha = 0.4) + 
  geom_line(data = quantmat_cond1, aes(y = p0.5), color = b_dark, size = 1) + 
  scale_x_continuous(breaks=pretty_breaks(10)) + 
  labs(x = "error (degrees) [red = pre, blue = post]", 
       y = "count +/- quantile", 
       subtitle = glue("per-condition posterior pred dist (median line, 90% interval over {nrow(post_predmean_sim)} sim datasets)\n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean draw)")) + 
  theme_cowplot()

#######################################################
# plot ecdf per condition

ggplot() +
  geom_ribbon(data = stim0_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "red", alpha = 0.3) + 
  geom_line(data = stim0_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "red", size = 1) +
  geom_ribbon(data = stim1_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "blue", alpha = 0.3) + 
  geom_line(data = stim1_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "blue", size = 1) + 
  scale_x_continuous(breaks=pretty_breaks(10)) + 
  geom_hline(yintercept = seq(0, 1, 0.25), linetype = "dashed", alpha = 0.2) +
  labs(x = "error (degrees) [red = pre, blue = post]", 
       y = "cumulative prob.", 
       subtitle = glue("per-condition posterior pred cdf (median line, 90% interval over {nrow(post_predmean_sim)} sim datasets) \n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean draw")) + 
  theme_cowplot()

* Notes

Predicting data, including group-level SD, shows that there is strong overlap in the pre and post predicitve dists.

Predicting data, ignoring group-level SD and predicting only the average subject, shows more of the trend of an effect.